Recurrence Quantification Analysis of Surface Electromyograms

In this paper we describe one of the possible application of recurrence plot strategy to the analysis of surface elctromyograms (sEMG). Recurrence quantification analysis (RQA) is an efficient time-series analysis tool pertaining to the class of non-linear dynamics time-domain processing. We analysed sEMG recorded during isometric contractions at constant and linearly varying force. In order to have a comparison with better known techniques, sEMG was also processed by standard FFT algorithms in order to compute the median frequency of sEMG power spectrum over time. From our results, we may conclude that RQA is almost equivalent to frequency domain analysis when considering costant force isometric contraction. At the contrary, it appears that RQA is most effective in detecting sEMG changes determined by brisk transients of force output. INTRODUCTION The analysis of the surface electromyographic signal (sEMG) is particularly attractive because it provides a relatively easy access to those physiological processes that allow the muscle to generate force and movement[1]. Up to now, the attempts to extract such information mainly relied on various approaches: • frequency domain analysis; • time domain analysis; • time/frequency analysis. The sEMG time course transformation into the frequency domain (sEMG) is a consolidated technique [2,3]. However, some assumptions on the characteristic of the sEMG and/or the experimental protocol must be made in order to apply the spectral analysis. In particular, a convenient epoch of signal (usually up to 1 second) must be taken into account in order to consider the signal as stationary; secondly, only isometric contractions, with few exceptions, are allowable to be studied. This latter point is based upon the necessity to minimise signal distortion introduced by relative movements between electrode and active muscle. Considering that sEMG may be represented as the product of a summation or superposition of the action potential trains of concurrently active motor units (MU), De Luca et al. [1] showed that the power density spectrum of the sEMG is obtained by summing all the auto and cross-spectra of the individual MUs action potential trains. From this representation these authors showed that the power spectrum contains information related to: [a] the shape of the motor units; [b] the mean firing rates of the individual MUs and [c] the time composition of the individual MU action potentials. In order to make easier the interpretation of spectrum modifications over time, the use of synthetic parameters obtained from the power spectrum has been proposed. In the literature (see [13] as an example), the most commonly used parameters are the mean frequency (MNF), the median frequency (MDF), the L/H ratio, the peak frequency, etc. All these parameters provide the common indication of a spectral compression toward the lower frequencies as the time of contraction increases. Among them, MNF and MDF are considered as good indicators of muscle fatigue. When the isometric force is not kept constant, but the subject is asked to vary the force during the motor task, the signal epoch width is suitably reduced in order to maintain the validity of the stationarity hypothesis. In this case, the time evolution of a given frequency parameter is considered to provide information related to MUs recruitment or de-recruitment. Concerning the time-domain analysis, one of the authors [4, 5] already studied the possibility to borrow from telecommunication science an algorithm for adaptive estimation able to extract dynamically the muscle force content from sEMG Usually, in order to describe various sEMG features during a given epoch, time domain parameters such as the root mean square (RMS), average rectified value (ARV), number of zero crossing are commonly used [1]. RMS and ARV provide an information which is affected by almost the same factors that affect the spectral parameters cited above. On the other hand, they are much sensitive to experimental errors (e.g. inter-electrode distance, electrode positioning, etc.) than spectral parameters [6]. By coupling the information obtained in the time and frequency domain during linearly varying isometric contraction, Bernardi et al [7] were able to show that, when recruitment is completed, the only way the system has to increase force is to increase the firing rate of the active MUs. Given the interest to study the sEMG in both time and frequency domains, recently some papers proposed different approaches using wavelets [8-12]. These are tools typically employed for non-stationary signals and are used to map them into the time/scale domain. Though the typical biomedical monodimensional wavelet-transform application is to the analysis of electroencephalograms (including the evoked responses) and electrocardiograms [13], in our opinion, the WT represents a very attractive approach to sEMG processing which will need further efforts. In the present paper, our intent is to describe the recurrence plot strategy which is one of the most efficient time-series analysis tools pertaining to the class of non-linear dynamics time-domain processing. Unlike other predominant time series techniques, recurrence plot analyses are not limited by data stationarity and/or size constraints. This graphical method, first introduced by Echmann et al in 1987 [14], was originally designed to locate hidden rhythms (the so-called recurring patterns) and nonstationarities (drifts) in experimental data sets. In particular, we illustrate how recurrence plots can take the sEMG signals, project them into multidimensional space by embedding procedures, and identify time recurrences (i.e., correlations) that are not so apparent in the original recordings. In accordance with Webber [15], the original description of recurrence plots is also extended by evaluating some specific recurrence variables (e.g., the percentage of determinism (%DET), the percentage of recurrence (%REC), etc) that quantify the deterministic structure and complexity of the plot itself. MATERIAL AND METHODS I) Subjects. Eleven voluntary healthy subjects (6 males and 5 females) gave their informed written consent to the research. At the time of experiments none was under regular training practice with the upper arm and all were right handed. Eight participated to isometric constant force (CF) experiment and three to isometric linearly varying (LV) force experiment. Electromyographyc signals were obtained from the right biceps brachii on both CF and LV experiments. II) Experimental set-up. Each subject seated comfortably having his right arm lying on a flat surface and his forearm in an hand made isometric device. The device was designed in order to precisely regulate the angle at elbow at 90°. A wrist strap was connected to a force transducer (Kistler, Type 9311A; Sensitivity = -3.93 pC/N, linearity < 0.3%) connected with a charge amplifier KISTLER 5001 (cut-off frequency at 180 kHz)). sEMG surface electrodes (3M) were Ag/AgCl 10 mm diameter on self adhesive support. The electrodes were positioned on the middle portion of muscle belly (short head) parallel to the longitudinal axis of muscle fibres, with an inter electrode distance (centre to centre) of 20 mm. sEMG signal was handled by GRASS amplifier with the pass band set between 1 Hz and 1 kHz. The acquisition system consisted of a multichannel PC board with a 12-bit successive approximation A/D converter (National Instruments, ATMIO-16E10). The sampling frequency was 2048 Hz. Calibration of the A/D converter and data acquisition were software controlled. Force data were also send through another AT-MIO-16E10 to a PC screen to be presented as target of the motor task to the subject. In this case sampling frequency was set to 10 samples per second. III)-Experimental protocol The biceps brachii (BB) muscle was studied in two different experimental situation: 1. Experiment A) -Constant force (CF) isometric contraction; 2. Experiment B)Linearly varying force (LV) isometric contraction. CF Experiment: Prior to the start of the actual experiment, the maximum voluntary contraction (MVC) was measured and the best of three attempts was taken as the MVC value. Then the subject was asked to produce three sets of contractions at 20%, 60% and 100% MVC. Each set was made by three attempts, each lasting 12 seconds. Each contraction was spaced by 10 minutes from the next in order to avoid fatigue phenomena. The recovery time was increased any time the subject was unable to match the 100% MVC. LV Experiment. This test was performed in order to study the effects of a non constant force isometric contraction on sEMG. During this task, which obviously introduced transients in the myoelectric signal, a computer generated a trapezoidal trajectory to be followed by the subject. The trajectory was designed in the following way. Three seconds at 0 force (preparatory phase); from 0 to 20%MVC in two seconds (slope 10% s-1); 15 seconds at 20% MVC; form 20% to 100%MVC in 8 seconds (slope 10% s-1); back to 20%MVC with the same slope. By applying the force to the transducer, the subject moved a crosshair cursor presented on the pC screen and tried to match at his best the pre-plotted trajectory. Force and sEMG data were sampled as in CF experiment. For any given experiment the recorded single differential signal was divided in epoch of 1 sec. In each epoch, the RQA and the Fourier analysis were applied and then the %DET and MDF were calculated and represented as a function of time. For CF experiment, absolute intercepts (INT) and slopes (SLO) of MDF and %DET were computed from linear regression over the entire duration of the contraction. . For LV experiments a 10th order polynomial interpolation was applied to the experimental points of MDF and %DET and the time instant in which the polynomial leaves out the confidence interval was detected. According to Webber et al [16], this time was considered the detected time of transient in the force level. II)-Recurrence quantification analysis (RQA) As stated above, the RQA method is able to detect state changes in drifting dynamical systems without necessitating any a-priori constraining assumptions on data set size, stationarity, or statistical distribution,etc. For these raisons, RQA is suitable for analysing sEMG during both constant force and non-constant force isometric contractions. While for the mathematical background we refer to the literature [14-17], here we point out the main choices within our application. Computations were performed on sEMG time-series data within episodic window consisting of M=2048 consecutive points corresponding to an epoch of 1 sec. Firstly, discrete data were lagged (by a number λ of samples equal to the first zero in the autocorrelation function) and embedded (Embedded Dimension ED=15), and then the absolute distance (DIST) between all possible vector pairs was computed to get the so called Embedding Matrix EM(i,j). While plotting the values of EM(i,j), for all values of i and j, would generate a three dimensional recurrence plot, viceversa plotting points at (i,j) co-ordinates for DIST values falling within a cut-off region (in our case a circular region with radius R=2) generates a two-dimensional recurrence plot. In Fig.1, where it is depicted the case of a sinusoid, a sinuosi+noise and a white noise recurrence plot, we may observe the ability of these technique to map the M size time-series into a bidimensional plot. In fact, the periodical characteristic of the sinusoid is in a way enhanced by means of that typical squared texture. If for periodic signals this feature is almost evident, at the contrary recurrence plots often contain subtle patterns that are not easily ascertained by qualitative visual inspection. For this reason, some quantitative descriptors able to emphasise different features of the plot have been introduced [16]. The main five parameteres are: percent determinism (%DET), percent recurrence